Code
source("../dsan-globals/_globals.r")PPOL 6805 / DSAN 6750: GIS for Spatial Data Science
Fall 2025
source("../dsan-globals/_globals.r")The Earth’s “width” is slightly greater than its “length” 😰






Using rnaturalearth with mapview
set.seed(6805)
library(tidyverse)── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.5
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ lubridate 1.9.4 ✔ tibble 3.3.0
✔ purrr 1.0.4 ✔ tidyr 1.3.1
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(sf)Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(rnaturalearth)
library(mapview)
france_sf <- ne_countries(country = "France", scale = 50)
(france_map <- mapview(france_sf, label = "geounit", legend = FALSE))france_cent_sf <- sf::st_centroid(france_sf)Warning: st_centroid assumes attributes are constant over geometries
france_map + mapview(france_cent_sf, label = "Centroid", legend = FALSE)Computing the union of all geometries in the sf via sf::st_union()
library(leaflet.extras2)Loading required package: leaflet
africa_sf <- ne_countries(continent = "Africa", scale = 50)
africa_union_sf <- sf::st_union(africa_sf)
africa_map <- mapview(africa_sf, label="geounit", legend=FALSE)
africa_union_map <- mapview(africa_union_sf, label="st_union(africa)", legend=FALSE)
africa_map | africa_union_mapafrica_bbox_sf <- sf::st_bbox(africa_sf)
africa_bbox_map <- mapview(africa_bbox_sf, label="st_bbox(africa)", legend=FALSE)
africa_map | africa_bbox_mapafrica_countries_cvx <- sf::st_convex_hull(africa_sf)
africa_countries_cvx_map <- mapview(africa_countries_cvx, label="geounit", legend=FALSE)
africa_map | africa_countries_cvx_mapUse st_union() first:
africa_cvx <- africa_sf |> st_union() |> st_convex_hull()
africa_cvx_map <- mapview(africa_cvx, label="geounit", legend=FALSE)
africa_map | africa_cvx_mapComputing the centroid of all geometries in the sf via sf::st_centroid()
africa_cents_sf <- sf::st_centroid(africa_sf)Warning: st_centroid assumes attributes are constant over geometries
africa_cents_map <- mapview(africa_cents_sf, label="geounit", legend=FALSE)
africa_map | africa_cents_mapnc <- system.file("shape/nc.shp", package="sf") |>
read_sf() |>
st_transform('EPSG:2264')
gr <- st_sf(
label = apply(expand.grid(1:10, LETTERS[10:1])[,2:1], 1, paste0, collapse = ""),
geom = st_make_grid(nc))
gr$col <- sf.colors(10, categorical = TRUE, alpha = .3)
# cut, to verify that NA's work out:
gr <- gr[-(1:30),]
suppressWarnings(nc_j <- st_join(nc, gr, largest = TRUE))
par(mfrow = c(2,1), mar = rep(0,4))
plot(st_geometry(nc_j), border = 'grey')
plot(st_geometry(gr), add = TRUE, col = gr$col)
text(st_coordinates(st_centroid(st_geometry(gr))), labels = gr$label, cex = .85)
# the joined dataset:
plot(st_geometry(nc_j), border = 'grey', col = nc_j$col)
text(st_coordinates(st_centroid(st_geometry(nc_j))), labels = nc_j$label, cex = .7)
plot(st_geometry(gr), border = '#88ff88aa', add = TRUE)
# Sample random points
africa_points_list <- sf::st_sample(africa_union_sf, 10)
africa_points_sf <- sf::st_sf(africa_points_list)
africa_points_map <- mapview(africa_points_sf, label="Random Point", col.regions=cb_palette[1], legend=FALSE)Warning in cbind(`Feature ID` = fid, mat): number of rows of result is not a
multiple of vector length (arg 1)
africa_map + africa_points_mapst_intersectscountries_w_points <- africa_sf[africa_points_sf,]
mapview(countries_w_points, label="geounit", legend=FALSE) + africa_points_maplengths()country_inter <- sf::st_intersects(africa_sf, africa_points_sf)
# Computes point counts for each polygon
(num_intersections <- lengths(country_inter)) [1] 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 0 0 0
[39] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 2
africa_sf <- africa_sf |> mutate(
num_points = num_intersections
) |> arrange(geounit)
africa_sf |> select(geounit, num_points) |> head()| geounit | num_points | geometry |
|---|---|---|
| Algeria | 2 | MULTIPOLYGON (((8.576563 36… |
| Angola | 2 | MULTIPOLYGON (((13.07275 -4… |
| Benin | 0 | MULTIPOLYGON (((1.622656 6…. |
| Botswana | 0 | MULTIPOLYGON (((25.25879 -1… |
| Burkina Faso | 0 | MULTIPOLYGON (((0.9004883 1… |
| Burundi | 0 | MULTIPOLYGON (((30.55361 -2… |
mapviewmapview(africa_sf, zcol="num_points")ggplot2Since we’re starting to get into data attributes rather than geometric features, switching to ggplot2 is recommended!
africa_sf |> ggplot(aes(fill=num_points)) +
geom_sf() +
theme_classic()
st_buffer()POINT) from Manhattan (POLYGON)?”POLYGON) / stretch of road (LINESTRING)?”POLYGONsKey line: manhattan_buffer_sf <- manhattan_union_sf |> st_buffer(dist = 1609.34) (1 mile \(\approx\) 1609.34 meters)
library(tidyverse)
library(sf)
library(tidycensus)
library(tigris)
library(mapview)
options(tigris_use_cache = TRUE)
manhattan_sf <- get_acs(
geography = "tract",
variables = "B19013_001",
state = "NY",
county = "New York",
year = 2020,
geometry = TRUE,
cb = FALSE
)
# Erase the island tracts real quick
island_tracts <- c(
"Census Tract 1, New York County, New York",
"Census Tract 2.02, New York County, New York"
)
manhattan_sf <- manhattan_sf |> filter(
!(NAME %in% island_tracts)
)
# Union of all census tracts within the county
manhattan_union_sf <- st_union(manhattan_sf)
manhattan_union_map <- mapview(manhattan_union_sf, label="New York County")
# Construct buffer (1 mile ~= 1609.34 meters)
manhattan_buffer_sf <- manhattan_union_sf |> st_buffer(dist = 1609.34)
manhattan_buffer_map <- mapview(manhattan_buffer_sf, label="Buffer (1 Mile)", col.regions = cbPalette[1], legend = TRUE)
manhattan_buffer_map + manhattan_union_mapLINESTRINGslibrary(tidyverse)
library(sf)
## st_buffer, style options (taken from rgeos gBuffer)
l1 = st_as_sfc("LINESTRING(0 0,1 5,4 5,5 2,8 2,9 4,4 6.5)")
op = par(mfrow=c(2,3))
plot(st_buffer(l1, dist = 1, endCapStyle="ROUND"), reset = FALSE, main = "endCapStyle: ROUND")
plot(l1,col='blue',add=TRUE)
plot(st_buffer(l1, dist = 1, endCapStyle="FLAT"), reset = FALSE, main = "endCapStyle: FLAT")
plot(l1,col='blue',add=TRUE)
plot(st_buffer(l1, dist = 1, endCapStyle="SQUARE"), reset = FALSE, main = "endCapStyle: SQUARE")
plot(l1,col='blue',add=TRUE)
plot(st_buffer(l1, dist = 1, nQuadSegs=1), reset = FALSE, main = "nQuadSegs: 1")
plot(l1,col='blue',add=TRUE)
plot(st_buffer(l1, dist = 1, nQuadSegs=2), reset = FALSE, main = "nQuadSegs: 2")
plot(l1,col='blue',add=TRUE)
plot(st_buffer(l1, dist = 1, nQuadSegs= 5), reset = FALSE, main = "nQuadSegs: 5")
plot(l1,col='blue',add=TRUE)
sf Documentationst_centroid()
POLYGON \(\mapsto\) POINTMULTIPOLYGON \(\mapsto\) POINTst_convex_hull()
POLYGON \(\mapsto\) POLYGONMULTIPOINT \(\mapsto\) POLYGONst_intersection()
POINT, POINT) \(\mapsto\) POINT | POINT EMPTYPOLYGON, POLYGON) \(\mapsto\) POLYGON | LINESTRING | POINT | POLYGON EMPTYst_is_empty() and st_dimension() become your new best friends 😉st_is_empty(): Distinguishes between, e.g., POINT EMPTY and POINT(0 0)st_dimension(): NA for empty versions, otherwise
2 for surfaces (POLYGON, MULTIPOLYGON)1 for lines (LINESTRING, MULTILINESTRING)0 for points (POINT, MULTIPOINT)
Unary Operations

Binary Operations

Ternary Operations

Quaternary Operations


(i spent 4 yrs of undergrad studying abstract algebra and now it all sits gathering dust somewhere deep within my brain plz just let me have this moment i’ll never mention it again i promise)

Each cell here visualizes one component of the DE-9IM string 1020F1102, which describes the relationship between the following two geometries:
POLYGON(0 0, 1 0, 1 1, 0 1, 0 0)LINESTRING(0.5 -0.5, 0.5 0.5)library(sf)
polygon <- po <- st_polygon(list(rbind(c(0,0), c(1,0), c(1,1), c(0,1), c(0,0))))
p0 <- st_polygon(list(rbind(c(-1,-1), c(2,-1), c(2,2), c(-1,2), c(-1,-1))))
line <- li <- st_linestring(rbind(c(.5, -.5), c(.5, 0.5)))
s <- st_sfc(po, li)
par(mfrow = c(3,3))
par(mar = c(1,1,1,1))
# "1020F1102"
# 1: 1
plot(s, col = c(NA, 'darkgreen'), border = 'blue', main = expression(paste("I(pol)",intersect(),"I(line) = 1")))
lines(rbind(c(.5,0), c(.5,.495)), col = 'red', lwd = 2)
points(0.5, 0.5, pch = 1)
# 2: 0
plot(s, col = c(NA, 'darkgreen'), border = 'blue', main = expression(paste("I(pol)",intersect(),"B(line) = 0")))
points(0.5, 0.5, col = 'red', pch = 16)
# 3: 2
plot(s, col = c(NA, 'darkgreen'), border = 'blue', main = expression(paste("I(pol)",intersect(),"E(line) = 2")))
plot(po, col = '#ff8888', add = TRUE)
plot(s, col = c(NA, 'darkgreen'), border = 'blue', add = TRUE)
# 4: 0
plot(s, col = c(NA, 'darkgreen'), border = 'blue', main = expression(paste("B(pol)",intersect(),"I(line) = 0")))
points(.5, 0, col = 'red', pch = 16)
# 5: F
plot(s, col = c(NA, 'darkgreen'), border = 'blue', main = expression(paste("B(pol)",intersect(),"B(line) = F")))
# 6: 1
plot(s, col = c(NA, 'darkgreen'), border = 'blue', main = expression(paste("B(pol)",intersect(),"E(line) = 1")))
plot(po, border = 'red', col = NA, add = TRUE, lwd = 2)
# 7: 1
plot(s, col = c(NA, 'darkgreen'), border = 'blue', main = expression(paste("E(pol)",intersect(),"I(line) = 1")))
lines(rbind(c(.5, -.5), c(.5, 0)), col = 'red', lwd = 2)
# 8: 0
plot(s, col = c(NA, 'darkgreen'), border = 'blue', main = expression(paste("E(pol)",intersect(),"B(line) = 0")))
points(.5, -.5, col = 'red', pch = 16)
# 9: 2
plot(s, col = c(NA, 'darkgreen'), border = 'blue', main = expression(paste("E(pol)",intersect(),"E(line) = 2")))
plot(p0 / po, col = '#ff8888', add = TRUE)
plot(s, col = c(NA, 'darkgreen'), border = 'blue', add = TRUE)
equals corresponds to the DE-9IM string "T*F**FFF*". If any two geometries obey this relationship, they are (topologically) equal!| 9IM | Interior | Boundary | Exterior |
|---|---|---|---|
| Interior | ![]() |
![]() |
![]() |
| Boundary | ![]() |
![]() |
![]() |
| Exterior | ![]() |
![]() |
![]() |
| 9IM | Interior | Boundary | Exterior |
|---|---|---|---|
| Interior | ![]() 2 |
![]() 1 |
![]() 2 |
| Boundary | ![]() 1 |
![]() 0 |
![]() 1 |
| Exterior | ![]() 2 |
![]() 1 |
![]() 2 |
| DE-9IM | Interior | Boundary | Exterior |
|---|---|---|---|
| Interior | 2 | 1 | 2 |
| Boundary | 1 | 0 | 1 |
| Exterior | 2 | 1 | 2 |
212101212

st_overlaps() |
Interior | Boundary | Exterior |
|---|---|---|---|
| Interior | T |
* |
T |
| Boundary | * |
* |
* |
| Exterior | T |
* |
* |
0, 1, 2):
T: “True” (non-empty, st_dimension() >= 0)F: “False” (empty, st_dimension() == NA)*: “Wildcard” (Don’t care what the value is)st_overlaps(): T*T***T**, st_equals(): T*F**FFF*| English | Mask | 212101212 |
Result |
|---|---|---|---|
| “Disjoint” | FF*FF**** |
FALSE |
x not disjoint from y |
| “Touches” | FT******* |
FALSE |
x doesn’t touch y |
| “Touches” | F***T**** |
FALSE |
x doesn’t touch y |
| “Crosses” | T*T***T** |
TRUE |
x crosses y |
| “Within” | TF*F***** |
FALSE |
x is not within y |
| “Overlaps” | T*T***T** |
TRUE |
x overlaps y |
st_relate(): The Ultimate Predicatelibrary(tidyverse)
library(rnaturalearth)
us <- ne_states(country="United States of America")
dc <- us |> filter(iso_3166_2 == "US-DC")
us <- us |>
mutate(
de9im = st_relate(us, dc),
touch = st_touches(us, dc, sparse = F)
) |>
select(iso_3166_2, name, de9im, touch) |>
arrange(name)although coordinates are longitude/latitude, st_relate assumes that they are
planar
us| iso_3166_2 | name | de9im | touch | geometry |
|---|---|---|---|---|
| US-AL | Alabama | FF2FF1212 | FALSE | MULTIPOLYGON (((-87.41958 3… |
| US-AK | Alaska | FF2FF1212 | FALSE | MULTIPOLYGON (((-141.0056 6… |
| US-AZ | Arizona | FF2FF1212 | FALSE | MULTIPOLYGON (((-111.0063 3… |
| US-AR | Arkansas | FF2FF1212 | FALSE | MULTIPOLYGON (((-90.30422 3… |
| US-CA | California | FF2FF1212 | FALSE | MULTIPOLYGON (((-114.7243 3… |
| US-CO | Colorado | FF2FF1212 | FALSE | MULTIPOLYGON (((-109.0463 4… |
| US-CT | Connecticut | FF2FF1212 | FALSE | MULTIPOLYGON (((-73.6417 41… |
| US-DE | Delaware | FF2FF1212 | FALSE | MULTIPOLYGON (((-75.05809 3… |
| US-DC | District of Columbia | 2FFF1FFF2 | FALSE | MULTIPOLYGON (((-77.02293 3… |
| US-FL | Florida | FF2FF1212 | FALSE | MULTIPOLYGON (((-87.44734 3… |
| US-GA | Georgia | FF2FF1212 | FALSE | MULTIPOLYGON (((-80.89029 3… |
| US-HI | Hawaii | FF2FF1212 | FALSE | MULTIPOLYGON (((-154.8996 1… |
| US-ID | Idaho | FF2FF1212 | FALSE | MULTIPOLYGON (((-117.0382 4… |
| US-IL | Illinois | FF2FF1212 | FALSE | MULTIPOLYGON (((-89.1237 36… |
| US-IN | Indiana | FF2FF1212 | FALSE | MULTIPOLYGON (((-84.80608 4… |
| US-IA | Iowa | FF2FF1212 | FALSE | MULTIPOLYGON (((-96.48266 4… |
| US-KS | Kansas | FF2FF1212 | FALSE | MULTIPOLYGON (((-102.0396 3… |
| US-KY | Kentucky | FF2FF1212 | FALSE | MULTIPOLYGON (((-89.42446 3… |
| US-LA | Louisiana | FF2FF1212 | FALSE | MULTIPOLYGON (((-89.52599 3… |
| US-ME | Maine | FF2FF1212 | FALSE | MULTIPOLYGON (((-71.08495 4… |
| US-MD | Maryland | FF2F11212 | TRUE | MULTIPOLYGON (((-75.64786 3… |
| US-MA | Massachusetts | FF2FF1212 | FALSE | MULTIPOLYGON (((-71.19396 4… |
| US-MI | Michigan | FF2FF1212 | FALSE | MULTIPOLYGON (((-84.4913 46… |
| US-MN | Minnesota | FF2FF1212 | FALSE | MULTIPOLYGON (((-97.22609 4… |
| US-MS | Mississippi | FF2FF1212 | FALSE | MULTIPOLYGON (((-88.40221 3… |
| US-MO | Missouri | FF2FF1212 | FALSE | MULTIPOLYGON (((-95.31725 4… |
| US-MT | Montana | FF2FF1212 | FALSE | MULTIPOLYGON (((-116.0482 4… |
| US-NE | Nebraska | FF2FF1212 | FALSE | MULTIPOLYGON (((-104.0537 4… |
| US-NV | Nevada | FF2FF1212 | FALSE | MULTIPOLYGON (((-114.0425 4… |
| US-NH | New Hampshire | FF2FF1212 | FALSE | MULTIPOLYGON (((-71.50585 4… |
| US-NJ | New Jersey | FF2FF1212 | FALSE | MULTIPOLYGON (((-75.54133 3… |
| US-NM | New Mexico | FF2FF1212 | FALSE | MULTIPOLYGON (((-108.1375 3… |
| US-NY | New York | FF2FF1212 | FALSE | MULTIPOLYGON (((-79.06523 4… |
| US-NC | North Carolina | FF2FF1212 | FALSE | MULTIPOLYGON (((-76.03173 3… |
| US-ND | North Dakota | FF2FF1212 | FALSE | MULTIPOLYGON (((-104.0476 4… |
| US-OH | Ohio | FF2FF1212 | FALSE | MULTIPOLYGON (((-80.52023 4… |
| US-OK | Oklahoma | FF2FF1212 | FALSE | MULTIPOLYGON (((-103.0002 3… |
| US-OR | Oregon | FF2FF1212 | FALSE | MULTIPOLYGON (((-124.4924 4… |
| US-PA | Pennsylvania | FF2FF1212 | FALSE | MULTIPOLYGON (((-79.76301 4… |
| US-RI | Rhode Island | FF2FF1212 | FALSE | MULTIPOLYGON (((-71.23686 4… |
| US-SC | South Carolina | FF2FF1212 | FALSE | MULTIPOLYGON (((-78.57316 3… |
| US-SD | South Dakota | FF2FF1212 | FALSE | MULTIPOLYGON (((-104.0567 4… |
| US-TN | Tennessee | FF2FF1212 | FALSE | MULTIPOLYGON (((-90.30422 3… |
| US-TX | Texas | FF2FF1212 | FALSE | MULTIPOLYGON (((-103.3115 2… |
| US-UT | Utah | FF2FF1212 | FALSE | MULTIPOLYGON (((-111.0502 4… |
| US-VT | Vermont | FF2FF1212 | FALSE | MULTIPOLYGON (((-73.35134 4… |
| US-VA | Virginia | FF2F11212 | TRUE | MULTIPOLYGON (((-76.01325 3… |
| US-WA | Washington | FF2FF1212 | FALSE | MULTIPOLYGON (((-122.753 48… |
| US-WV | West Virginia | FF2FF1212 | FALSE | MULTIPOLYGON (((-82.58945 3… |
| US-WI | Wisconsin | FF2FF1212 | FALSE | MULTIPOLYGON (((-87.80425 4… |
| US-WY | Wyoming | FF2FF1212 | FALSE | MULTIPOLYGON (((-109.0463 4… |
us |> filter(touch)Warning: Using one column matrices in `filter()` was deprecated in dplyr 1.1.0.
ℹ Please use one dimensional logical vectors instead.
| iso_3166_2 | name | de9im | touch | geometry |
|---|---|---|---|---|
| US-MD | Maryland | FF2F11212 | TRUE | MULTIPOLYGON (((-75.64786 3… |
| US-VA | Virginia | FF2F11212 | TRUE | MULTIPOLYGON (((-76.01325 3… |
library(mapview)
library(leaflet.extras2)
africa_sf <- ne_countries(continent = "Africa", scale = 50) |> select(iso_a3, geounit)
africa_map <- mapview(africa_sf, label="geounit", legend=FALSE)
N <- 10
africa_union_sf <- sf::st_union(africa_sf)
sampled_points_sf <- sf::st_sample(africa_union_sf, N) |> sf::st_sf() |> mutate(temp = runif(N, 0, 100))
sampled_points_map <- mapview(sampled_points_sf, label="Random Point", col.regions=cbPalette[1], legend=FALSE)
countries_points_sf <- africa_sf[sampled_points_sf,]
filtered_map <- mapview(countries_points_sf, label="geounit", legend=FALSE) + sampled_points_map
(africa_map + sampled_points_map) | filtered_mapThe issue: Data attributes of POINTs are not merged into data attributes of POLYGONs
POINT Attributes
st_geometry(sampled_points_sf) <- c("geom")
sampled_points_sf |> head()| geom | temp |
|---|---|
| POINT (-5.135802 22.30373) | 99.14195 |
| POINT (48.01641 10.55753) | 41.21422 |
| POINT (47.30147 5.038318) | 95.30175 |
| POINT (22.64591 7.799242) | 68.37119 |
| POINT (7.739425 5.95728) | 56.40943 |
| POINT (22.13408 20.87522) | 86.12631 |
POLYGON Attributes
countries_points_sf |> head(4)| iso_a3 | geounit | geometry | |
|---|---|---|---|
| 2 | ZMB | Zambia | MULTIPOLYGON (((30.39609 -1… |
| 58 | SOM | Somalia | MULTIPOLYGON (((41.53271 -1… |
| 59 | -99 | Somaliland | MULTIPOLYGON (((48.93857 11… |
| 91 | NGA | Nigeria | MULTIPOLYGON (((7.300781 4…. |
st_join()joined_sf <- countries_points_sf |> st_join(sampled_points_sf)
joined_sf |> head()| iso_a3 | geounit | temp | geometry | |
|---|---|---|---|---|
| 2 | ZMB | Zambia | 5.788125 | MULTIPOLYGON (((30.39609 -1… |
| 58 | SOM | Somalia | 95.301751 | MULTIPOLYGON (((41.53271 -1… |
| 59 | -99 | Somaliland | 41.214222 | MULTIPOLYGON (((48.93857 11… |
| 91 | NGA | Nigeria | 56.409432 | MULTIPOLYGON (((7.300781 4…. |
| 114 | MLI | Mali | 99.141948 | MULTIPOLYGON (((-11.3894 12… |
| 123 | LBY | Libya | 86.126311 | MULTIPOLYGON (((9.51875 30…. |

g <- st_make_grid(st_bbox(st_as_sfc("LINESTRING(0 0,1 1)")), n = c(2,2))
par(mar = rep(0,4))
plot(g)
plot(g[1] * diag(c(3/4, 1)) + c(0.25, 0.125), add = TRUE, lty = 2)
text(c(.2, .8, .2, .8), c(.2, .2, .8, .8), c(1,2,4,8), col = 'red')
Assume the extensive attribute \(Y\) is uniformly distributed over a space \(S_i\) (e.g., for population counts we assume everyone is evenly-spaced across the region)
We first compute \(Y_{ij}\), derived from \(Y_i\) for a sub-area of \(S_i\), \(A_{ij} = S_i \cap T_j\):
\[ \hat{Y}_{ij}(A_{ij}) = \frac{|A_{ij}|}{|S_i|}Y_i(S_i) \]
where \(|\cdot|\) denotes area.
Then we can compute \(Y_j(T_j)\) by summing all the elements over area \(T_j\):
\[ \hat{Y}_j(T_j) = \sum_{i=1}^{p}\frac{|A_{ij}|}{|S_i|}Y_i(S_i) \]
\[ \hat{Y}_{ij} = Y_i(S_i) \]
\[ \hat{Y}_j(T_j) = \sum_{i=1}^{p}\frac{|A_{ij}|}{|T_j|}Y_j(S_i) \]